Iago Giné Vázquez
2021-10-26/27
>. Ahí se escriben las instrucciones. Se ejecutan con Enter..r en lugar de .txt..r en el que se puede escribir el código (con autocompletado) y ejecutarloControl+Enter"a", o 'z', hasta un texto largo como "Lorem ipsum dolor sit amet..." o 'Lorem ipsum dolor sit amet...'f seguido de los objetos o valores que necesita entre paréntesis.
f()x con una función f se codifica f(x).#getwd;setwd evaluando la dirección del directorio de trabajo en que queremos situarnos ahora;list.files;getwd() # el directorio donde R está trabajando
setwd("P:/projects/Curso-R") # (en Windows) se tienen que sustituir los \ por /
list.files()O bien, mediante RStudio:
`Session > Set Working Directory > Choose Directory`
`Session > Set Working Directory > To Source File Location`
tidyverse evaluamos con la función install.packages la cadena consistente en el nombre de la librería, "tidyverse";readxl, evaluamos con la función library (o bien require, una función distinta), bien la cadena con el nombre de la librería "readxl" o bien (el nombre de) la librería misma, readxl.install.packages("tidyverse") # instala la librería tidyverse, que a su vez instala las librerías readxl, dplyr, tidyr, haven, y otras (https://www.tidyverse.org, https://github.com/tidyverse/tidyverse)
library(readxl)
# library("readxl")
# require(readxl)
# require("readxl")
# (.packages()) # vemos qué librerías tenemos cargadas en este momentoO bien, mediante RStudio:
Panel `Packages > Install`
help() o mediante ?. Los siguientes ejemplos nos muestran los documentos de ayuda de las respectivas funciones.library(readxl) asignamos (el nombre de) la librería readxl al primer argumento de la función library, que es package. Entonces, otra forma de ejecutar library(readxl) es ejecutar library(package = readxl). En este caso, es exactamente lo mismo, pero cuando queremos usar otro argumento de una función distinto del primero en general tenemos que nombrarlo.?help
help(package = "readxl") # Informémonos acerca de la librería readxl
?read_excel # ayuda de la función read_excelPara abrir ficheros de excel disponemos de diversas funciones en las librerías readxl y openxlsx entre muchas otras. Aunque la segunda librería tiene más funciones y opciones para manipular y escribir ficheros .xlsx, la primera dispone de funciones que nos permiten leer también ficheros .xls. En cualquier caso, estas funciones evalúan, bien la dirección absoluta del fichero que hay que abrir, bien la dirección relativa con respecto al directorio de trabajo de R como en el siguiente ejemplo.
## # A tibble: 4,753 × 12
## q0002_hhid sex age mar_1 edu1 phys_hea1 hea1 dep1 score_lon1 score_sup1 income1 income_inf1
## <dbl> <chr> <chr> <chr> <chr> <chr> <dbl> <chr> <dbl> <dbl> <dbl> <chr>
## 1 1 fem 35-49 Yes Tertiary None 70.5 No 3 12 15 Normal
## 2 2 masc 35-49 Yes Tertiary None 78.9 No 3 12 17 Normal
## 3 3 masc 65-79 Yes Primary None 66.6 No 3 12 15 Normal
## 4 4 masc 50-64 Yes Primary 1 physical health problem 77.4 No 3 11 18 Normal
## 5 5 fem 50-64 No Tertiary None 52.5 No 3 12 29 Normal
## 6 6 fem 50-64 No Tertiary None 53.7 No 3 13 25 Normal
## 7 7 masc 50-64 Yes Secondary None 59.1 No 6 10 12 Good
## 8 8 masc 50-64 Yes Primary None 76.4 No 3 12 14 Normal
## 9 9 fem 50-64 No Secondary None 73.7 No 3 11 26 Good
## 10 10 fem 80+ No Primary 2+ physical health problems 28.0 No 6 13 10 Normal
## # … with 4,743 more rows
Nota: Para abrir ficheros *.csv disponemos de la función read.csv en R, de la función read_csv en la librería readr y de la función fread en la librería data.table entre otras. La librería haven dispone de diversas funciones que permiten abrir y escribir ficheros en formatos de Stata, SPSS y SAS (pueden consultarse en su ayuda, package(help = "haven")). Por otra parte la librería readspss (https://github.com/JanMarvin/readspss) tiene funciones que permiten abrir bases de datos en formato de SPSS encriptadas con contraseña. Podemos encontrar más información en esta breve introducción a R
La función read_excel abre la base de datos, pero a diferencia de instrucciones como use o import excel en Stata, no la copia a memoria. Así, si en Stata podríamos ejecutar describe o summarize, si en R ejecutamos una función como summary el resultado es la propia función.
## function (object, ...)
## UseMethod("summary")
## <bytecode: 0x561722765af8>
## <environment: namespace:base>
Tenemos que evaluar la base de datos que abrimos con read_excel con la función summary:
## q0002_hhid sex age mar_1 edu1 phys_hea1 hea1 dep1 score_lon1 score_sup1 income1 income_inf1
## Min. : 1 Length:4753 Length:4753 Length:4753 Length:4753 Length:4753 Min. : 0.00 Length:4753 Min. :3.00 Min. : 3.00 Min. : 1.000 Length:4753
## 1st Qu.: 1327 Class :character Class :character Class :character Class :character Class :character 1st Qu.:42.83 Class :character 1st Qu.:3.00 1st Qu.:10.00 1st Qu.: 2.000 Class :character
## Median : 2615 Mode :character Mode :character Mode :character Mode :character Mode :character Median :55.64 Mode :character Median :3.00 Median :12.00 Median : 5.000 Mode :character
## Mean : 3145 Mean :53.97 Mean :3.74 Mean :11.53 Mean : 8.989
## 3rd Qu.: 3807 3rd Qu.:66.06 3rd Qu.:4.00 3rd Qu.:13.00 3rd Qu.:15.000
## Max. :92755 Max. :92.82 Max. :9.00 Max. :14.00 Max. :35.000
## NA's :187 NA's :285 NA's :414 NA's :1067
Abrir la base de datos mediante la función read_excel cada vez que queremos operar con ella no es práctico. En lugar de esto, guardamos la base de datos en una variable.
R es un lenguaje de programación que nos permite acceder a objetos con distintas estructuras y manipularlos. Los tipos de objetos más importantes son:
TRUE, FALSE y NA (missing)clistNULLLa mayoría de los objetos pueden tener atributos, como pueden ser las dimensiones o las etiquetas. De esta manera se pueden crear objetos con estructuras más complejas:
dim incluyendo las 2 (resp. -) dimensiones del objeto. Podemos crear matrices mediante la función matrix.levels que corresponde a las posibles categorías de una variable categórica y un atributo class="factor". Por defecto son nominales, pero también pueden tener categorías ordenadas, caso en el que el atributo class es el vector c("ordered"," factor")?DateTimeClasses puede verse más información)
POSIXct (resp. Date)POSIXltnames que designa los nombres de las columnas o variables de la base de datos y un atributo row.names que designa las filas u observaciones. Podemos crearlos mediante la función data.frameDisponemos de un conjunto de operadores en R que nos permiten acceder y manipular los objetos. Además de los operadores aritméticos y lógicos, entre otros disponemos de los siguientes:
<-, para asignar un nombre a un objeto. Por ejemplo, x <- c(3,log(4), 5.2), y <- list(3, 'a', list(TRUE)), M <- matrix(c(1,2,3,4), nrow = 2, ncol = 2, byrow=TRUE)[, para acceder a un elemento de un vector o a una sublista. Por ejemplo, x[3], y[2], y[3], M[,2] o M[1,2][[, para acceder a un elemento en una lista. Por ejemplo, y[[2]]NA, existente para los distintos tipos básicos, o también NaN en el caso de objetos numéricos.nombre = elemento. Por ejemplo y <- list(numero = 3, `letra` = 'a', "sublista" = list(TRUE)).y[c("letra", 'sublista')].[[ o $. Por ejemplo y[['letra']] o y[["sublista"]], o equivalentemente y$letra o y$sublista (o y$`letra` o y$`sublista`).print.length, class, str y attributes respectivamente.names nos permite conocer los nombres de los elementos de un objeto, como pueden ser las variables de un data frame.levels evaluada en un factor retorna las categorías del factor.dim vemos las dimensiones de un objeto, como pueden ser las filas y columnas de una matriz o de un data frame. Para objetos similares, dimnames informa acerca de los nombres de las filas y columnas.is y as.
is permiten verificar que un objeto es de una determinada clase. Algunos ejemplos:
is.na(y)is.character(y$numero), is.numeric(y$letra)is.data.frame(y)as fuerzan a los objetos a transformarse a la clase especificada cuando esto es posible.
as.character(y)as.numeric(list("3", TRUE, FALSE, "a"))as.Date("01/05/1965", "%m/%d/%Y")as.data.frame(M)seq genera una sucesión de números desde uno inicial hasta uno final, y permite especificar o el tamaño de los intervalos o la longitud de la sucesión. Por ejemplo seq(from = 5, to = 20, by = 2): también permite generar sucesiones de números enteros. Por ejemplo, 5:20rep genera un vector de elementos a partir de un vector inicial dado. Por ejemplo, rep(x = c(1,3,5), times = 2) o rep(c(1,3,5), each = 5).paste y paste0 combinan cadenas de caracteres según se especifique. Por ejemplo, paste(y$letra, c(1,2))Asignamos la base de datos a una variable que llamamos dataw1. A partir de este momento tenemos un objeto de R con nombre dataw1, ya no es una base de datos en excel.
dataw1 <- read_excel("data_curs_stat/EP1.xls", sheet = 1) # aunque en este caso no es necesario, especificamos también que abrimos la hoja 1 del fichero excel
class(dataw1)## [1] "tbl_df" "tbl" "data.frame"
## [1] 12
## [1] "q0002_hhid" "sex" "age" "mar_1" "edu1" "phys_hea1" "hea1" "dep1" "score_lon1" "score_sup1" "income1" "income_inf1"
## tibble [4,753 × 12] (S3: tbl_df/tbl/data.frame)
## $ q0002_hhid : num [1:4753] 1 2 3 4 5 6 7 8 9 10 ...
## $ sex : chr [1:4753] "fem" "masc" "masc" "masc" ...
## $ age : chr [1:4753] "35-49" "35-49" "65-79" "50-64" ...
## $ mar_1 : chr [1:4753] "Yes" "Yes" "Yes" "Yes" ...
## $ edu1 : chr [1:4753] "Tertiary" "Tertiary" "Primary" "Primary" ...
## $ phys_hea1 : chr [1:4753] "None" "None" "None" "1 physical health problem" ...
## $ hea1 : num [1:4753] 70.5 78.9 66.6 77.4 52.5 ...
## $ dep1 : chr [1:4753] "No" "No" "No" "No" ...
## $ score_lon1 : num [1:4753] 3 3 3 3 3 3 6 3 3 6 ...
## $ score_sup1 : num [1:4753] 12 12 12 11 12 13 10 12 11 13 ...
## $ income1 : num [1:4753] 15 17 15 18 29 25 12 14 26 10 ...
## $ income_inf1: chr [1:4753] "Normal" "Normal" "Normal" "Normal" ...
## # A tibble: 6 × 12
## q0002_hhid sex age mar_1 edu1 phys_hea1 hea1 dep1 score_lon1 score_sup1 income1 income_inf1
## <dbl> <chr> <chr> <chr> <chr> <chr> <dbl> <chr> <dbl> <dbl> <dbl> <chr>
## 1 1 fem 35-49 Yes Tertiary None 70.5 No 3 12 15 Normal
## 2 2 masc 35-49 Yes Tertiary None 78.9 No 3 12 17 Normal
## 3 3 masc 65-79 Yes Primary None 66.6 No 3 12 15 Normal
## 4 4 masc 50-64 Yes Primary 1 physical health problem 77.4 No 3 11 18 Normal
## 5 5 fem 50-64 No Tertiary None 52.5 No 3 12 29 Normal
## 6 6 fem 50-64 No Tertiary None 53.7 No 3 13 25 Normal
$ y [[## [1] "character"
#summary(dataw1$phys_hea1)
str(dataw1$phys_hea1) # chr quiere decir que phys_hea1 es una variable carácter, es decir, un vector de tipo cadena de carácteres## chr [1:4753] "None" "None" "None" "1 physical health problem" "None" "None" "None" "None" "None" "2+ physical health problems" "None" "None" "2+ physical health problems" "None" NA NA "2+ physical health problems" "None" ...
## NULL
##
## 1 physical health problem 2+ physical health problems None
## 1538 887 1723
##
## 1 physical health problem 2+ physical health problems None <NA>
## 1538 887 1723 605
##
## 1 physical health problem 2+ physical health problems None
## 0.3707811 0.2138380 0.4153809
Es posible introducir diferentes instrucciones en una misma línea separadas por ; (indiferentemente de los espacios en blanco en medio)
## [1] "numeric"
## num [1:4753] 70.5 78.9 66.6 77.4 52.5 ...
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 42.83 55.64 53.97 66.06 92.82 187
+, -, *, /, %% (residuo), ^ (potencia)sum, round, floor, log, log10, log2, exp, sqrt, factorial== (igualdad), ! (negación), != (no igualdad), <, <=, >, >= (desigualdades), %in% (pertenencia), |, || (o), &, && (y)## [1] TRUE
## [1] FALSE
## [1] FALSE
## [1] TRUE
## [1] TRUE
## [1] FALSE
## [1] FALSE
## [1] TRUE
## [1] 2
openxlsx y buscar la ayuda de sus funciones read.xlsx y loadWorkbookpaste, seq o rep) qué opciones presentan según como se especifiquen sus argumentos. Para qué sirven las funciones seq_len y seq_along?setwd y read_excel convenientemente (hay varias posibilidades), guardar también las bases de datos EP2.xls y EP3.xls en dos variables que llamaremos dataw2 y dataw3. Podéis ver dataw1, dataw2, dataw3 y sus dimensiones en el panel Environment de RStudio?c tiene todos los elementos del mismo tipo y que NA es un elemento de tipo lógico, por qué la instrucción c(1,3,5,7, NA, 8,5, NA) no da lugar a un error?| y || , y & y && se diferencian cuando comparan vectores lógicos con más de un elemento, como por ejemplo c(TRUE, FALSE). Cuáles son esas diferencias?Nota: para buscar la ayuda de un operador como | podemos ejecutar la instrucción ?`|`
as.factor permite transformar a un factor, objeto que codifica una variable categórica:## [1] "factor"
## Factor w/ 3 levels "1 physical health problem",..: 3 3 3 1 3 3 3 3 3 2 ...
## $levels
## [1] "1 physical health problem" "2+ physical health problems" "None"
##
## $class
## [1] "factor"
##
## 1 physical health problem 2+ physical health problems None <NA>
## 1538 887 1723 605
## [1] "1 physical health problem" "2+ physical health problems" "None"
Si nos interesa un atributo concreto, usamos la función attr:
## [1] "factor"
## [1] "1 physical health problem" "2+ physical health problems" "None"
Del mismo modo que podemos encontrarnos con variables categóricas ordenadas, como es el caso de las escalas de tipo Likert, los factores pueden ser ordenados. En este caso podemos utilizar la función as.ordered, aunque ésta puede ordenar las categorías erróneamente.
## $levels
## [1] "1 physical health problem" "2+ physical health problems" "None"
##
## $class
## [1] "ordered" "factor"
## Ord.factor w/ 3 levels "1 physical health problem"<..: 3 3 3 1 3 3 3 3 3 2 ...
Ante ello, es preferible construír el factor especificando el orden que nos interesa:
x <- factor(dataw1$phys_hea1, ordered = TRUE, levels = c("None", "1 physical health problem", "2+ physical health problems"))
str(x)## Ord.factor w/ 3 levels "None"<"1 physical health problem"<..: 1 1 1 2 1 1 1 1 1 3 ...
## $levels
## [1] "None" "1 physical health problem" "2+ physical health problems"
##
## $class
## [1] "ordered" "factor"
## [1] "None" "1 physical health problem" "2+ physical health problems"
## x
## None 1 physical health problem 2+ physical health problems <NA>
## 1723 1538 887 605
dplyr)Guardamos la transformación en la base de datos. Para ello usaremos la función mutate de la librería dplyr, que permite realizar varias transformaciones separadas por comas
library(dplyr)
#?mutate
dataw1 <- mutate(dataw1, phys_hea1 = as.factor(phys_hea1))
# es lo mismo que:
dataw1 <- dataw1 %>% mutate(phys_hea1 = as.factor(phys_hea1))
# podemos escrbir el argumento a la derecha de %>% en las líneas inferiores
dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
mutate(phys_hea1 = as.factor(phys_hea1)) # transformamos la variable y la guardamos con el mismo nombre
# finalmente la base de datos queda guardada con el mismo nombre mediante la asignación inicial (dataw1 <- ...)mutate es la base de datos o data frame.variable_output = transformación separadas por comas.%>%, la función situada en el lado derecho (en este caso mutate) ya reconoce que su primer argumento es el valor situado en el lado izquierdo de %>% (en este caso dataw1)dplyr porque incluye otras funciones que facilitan realizar diversas operaciones y que veremos más adelante, de manera que usemos un lenguaje consistente. Sin embargo, podríamos escribir la instrucción especificada encima sin acudir a ninguna librería:# Forma simple
dataw1$phys_hea1 <- as.factor(dataw1$phys_hea1)
# Usando la función transform
dataw1 <- transform(dataw1, phys_hea1 = as.factor(phys_hea1))Y a partir de la versión 4.1.0, R dispone del nuevo operador |>:
## [1] "q0002_hhid" "sex" "age" "mar_1" "edu1" "phys_hea1"
dataw1 %>% # la pipa más habitual: envía los datos de la izquierda a la función que sigue (a la derecha o abajo) y retorna su evaluación
names() %>% # esta pipa está incluida en la librería dplyr, pero es original de la librería magrittr
head()## [1] "q0002_hhid" "sex" "age" "mar_1" "edu1" "phys_hea1"
## phys_hea1
## edu1 1 physical health problem 2+ physical health problems None
## Less primary 445 396 304
## Primary 493 276 471
## Secondary 406 169 625
## Tertiary 192 46 320
Nota: la pipa %>% puede ser problemática con algunas funciones, como por ejemplo save. La pipa de R base |>, aunque no se puede utilizar para determinadas estructuras más complejas, es más segura. Por ejemplo, dataw1 |> save(file = "datos.Rdata") es exactamente lo mismo que save(dataw1, file = "datos.Rdata), pero dataw1 %>% save(file = "datos.Rdata") implica una instrucción ligeramente distinta.
dplyrComo la variable q0002_hhid es un id podría interesarnos que fuera de clase cadena
Varias transformaciones seguidas las podemos evaluar:
mutate:dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
mutate(phys_hea1 = as.factor(phys_hea1)) %>% # transformamos la variable phys_hea1 y la guardamos con el mismo nombre, y entonces
mutate(dep1 = as.factor(dep1)) %>% # transformamos la variable dep1 y la guardamos con el mismo nombre, y entonces
mutate(q0002_hhid = as.character(q0002_hhid)) # transformamos la variable q0002_hhid y la guardamos con el mismo nombremutate:dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
mutate(phys_hea1 = as.factor(phys_hea1), dep1 = as.factor(dep1), q0002_hhid = as.character(q0002_hhid)) # transformamos las variables phys_hea1, dep1 y q0002_hhid y las guardamos con los mismos nombresas.factor) a varias variables (y guardándolas con el mismo nombre), podemos usar across(vector de variables, función) dentro de mutate:dataw1 <- dataw1 %>% #partimos de la base de datos dataw1 y entonces
mutate(across(c(phys_hea1, dep1), as.factor), q0002_hhid = as.character(q0002_hhid)) # evaluamos as.factor a través de las variables phys_hea1 y dep1 y transformamos también number_idNota: Dentro de la función across, podemos seleccionar las variables que nos interesan de formas alternativas a escribir un vector de variables como c(phys_hea1, dep1). Veremos una introducción a estas alternativas en la siguiente diapositiva y en los ejercicios
dplyr: selección de variablesSelección de variables
dataw1 %>% #partimos de la base de datos dataw1 y entonces
select(q0002_hhid, hea1, dep1) %>% # mantenemos sólo las columnas q0002_hhid, hea1 y dep1, y entonces
names() # miramos qué nombres tienen## [1] "q0002_hhid" "hea1" "dep1"
dataw1 %>%
select(ends_with("1")) %>% # seleccionamos todas las variables cuyo nombre acaba en 1, y entonces
names() # miramos qué nombres tienen## [1] "mar_1" "edu1" "phys_hea1" "hea1" "dep1" "score_lon1" "score_sup1" "income1" "income_inf1"
Nota: Podemos seleccionar variables por índice y por nombre. Por ejemplo, para las 3 primeras variables escribiríamos select(1:3).
Pull. Una función en dplyr para el operador $
## Factor w/ 3 levels "1 physical health problem",..: 3 3 3 1 3 3 3 3 3 2 ...
Posibles combinaciones distintas de edad, educación y depresión
## # A tibble: 59 × 3
## age edu1 dep1
## <chr> <chr> <chr>
## 1 35-49 Tertiary No
## 2 65-79 Primary No
## 3 50-64 Primary No
## 4 50-64 Tertiary No
## 5 50-64 Secondary No
## 6 80+ Primary No
## 7 18-34 Tertiary No
## 8 18-34 Secondary No
## 9 65-79 Secondary <NA>
## 10 50-64 Tertiary <NA>
## # … with 49 more rows
dplyr y funciones de R para el cálculo de estadísticosdataw1 %>% #partimos de la base de datos dataw1 y entonces
select(q0002_hhid, hea1, dep1) %>% # mantenemos sólo las columnas q0002_hhid, hea1 y dep1, y entonces
filter(dep1 == "Yes") %>% # nos quedamos con las filas de quienes padecen depresión, y entonces
arrange(desc(hea1)) %>% # ordenamos las filas por de mayor a menor valor de estado de salud, y entonces
head(3) # nos quedamos con los 3 pacientes con depresión con mayor valor de Health state## # A tibble: 3 × 3
## q0002_hhid hea1 dep1
## <dbl> <dbl> <chr>
## 1 4941 82.1 Yes
## 2 4339 80.2 Yes
## 3 1299 79.8 Yes
mean, median o quantile; se pueden usar otras como var, sd, min, max, IQR, etc.)dataw1 %>% #partimos de la base de datos dataw1 y entonces
select(q0002_hhid, hea1, dep1, score_lon1) %>% # mantenemos sólo las columnas number_id, hea1 y dep1, y entonces
group_by(dep1) %>% # agrupamos por las categorías de depresión, y entonces
summarise(mean_hea1 = mean(hea1, na.rm = TRUE), median_hea1 = median(hea1, na.rm = TRUE), tercil_hea1 = quantile(hea1, probs = 1/3, na.rm = TRUE), n = n(), distinct_score_lon1 = n_distinct(score_lon1)) ## # A tibble: 3 × 6
## dep1 mean_hea1 median_hea1 tercil_hea1 n distinct_score_lon1
## <chr> <dbl> <dbl> <dbl> <int> <int>
## 1 No 55.9 57.4 50.2 4073 8
## 2 Yes 38.4 37.9 31.1 510 8
## 3 <NA> NaN NA NA 170 1
# para cada categoría de depresión calculamos media, mediana y primer tercil de hea1,
# número de observaciones y número de observaciones distintas de score_lon1
dataw1 %>% distinct(dep1, score_lon1) %>% group_by(dep1) %>% summarise(n = n())## # A tibble: 3 × 2
## dep1 n
## <chr> <int>
## 1 No 8
## 2 Yes 8
## 3 <NA> 1
Nota: usamos na.rm = TRUE dentro de mean, de median y de quantile para que calcule la media de aquellos valores que no son missing. En caso contrario, cuando hay missings el resultado es NA.
Nota: Todo lo anterior se puede realizar también con funciones de R sin necesidad de acudir a la librería dplyr, pero las alternativas pueden ser más complejas.
Prevalencia de depresión:
dataw1 %>%
count(dep1) %>% # Contamos los individuos en cada categoría de depresión y entonces
mutate(prop = 100*n/sum(n)) # calculamos el %## # A tibble: 3 × 3
## dep1 n prop
## <chr> <int> <dbl>
## 1 No 4073 85.7
## 2 Yes 510 10.7
## 3 <NA> 170 3.58
Prevalencia de depresión por grupos de edad y sin tener en cuenta los missings:
# ?is.na
dataw1 %>%
group_by(age) %>%
count(dep1) %>% # Contamos los individuos en cada categoría de depresión por grupo de edad y entonces
filter(!is.na(dep1)) %>% # eliminamos los missings si nos interesa contar el porcentaje sobre el total de respuestas válidas, y entonces
mutate(prop = 100*n/sum(n)) # Calculamos el %## # A tibble: 10 × 4
## # Groups: age [5]
## age dep1 n prop
## <chr> <chr> <int> <dbl>
## 1 18-34 No 382 93.9
## 2 18-34 Yes 25 6.14
## 3 35-49 No 500 90.7
## 4 35-49 Yes 51 9.26
## 5 50-64 No 1558 88.5
## 6 50-64 Yes 202 11.5
## 7 65-79 No 1298 87.3
## 8 65-79 Yes 188 12.7
## 9 80+ No 335 88.4
## 10 80+ Yes 44 11.6
Hemos dicho que R es como una calculadora, y que si no asignamos los objetos a variables, se muestran en consola pero no quedan guardados en ningún sitio.
Environment de RStudio o podemos ver sus nombres todavía con la función ls. Se pueden eliminar objetos mediante la función rm, separando sus nombres por comas.save en ficheros con la extensión .rda o .rdata (aunque a veces también se escribe la r e incluso la d en mayúscula, por ejemplo .RData).load (o de forma distinta, con la función attach)..csv o con formatos de excel. De hecho, a menudo las funciones mencionadas de la forma read_ o read. tienen funciones correspondientes write_ o write.. Por ejemplo, podríamos guardar el data frame dataw1 en un fichero csv con write.csv:sink. El primer argumento de esta función es el nombre del fichero donde se guardarán los resultados.getwd())#?sink
sink("Summary_Estudi_poblacional_w1.txt", split = TRUE) # El fichero de texto donde se guardarán los resultados se llamará Summary_Estudi_poblacional_w1.txt,
# pero podría llamarse como cualquier otro fichero de texto, no hay más restricción para el nombre que las del sistema operativo donde se está trabajando
table(dataw1$phys_hea1)##
## 1 physical health problem 2+ physical health problems None
## 1538 887 1723
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 42.83 55.64 53.97 66.06 92.82 187
q0002_hhid a cadena de carácteres.- delante de end_with en dataw1 %>% select(ends_with("1")) %>% names()?ends_with disponemos?dataw1 %>% distinct(dep1, score_lon1) %>% group_by(dep1) %>% summarise(n = n()) y las distintas funciones mencionadas referentes a estadísticos. Mirar también qué hace y cómo se puede usar la función rangeif, else, ifelse (e if_else, la función equivalente en dplyr), case_when (otra función en dplyr)for, while## [1] "Tertiary" "Tertiary" "Primary" "Primary" "Tertiary" "Tertiary" "Secondary" "Primary" "Secondary" "Primary"
## x
## Less primary Primary Secondary Tertiary <NA>
## 1364 1426 1329 629 5
for(i in seq_along(x)){
if(is.na(x[i])){
y[i] <- "Unknown"
} else if (x[i] == "Less primary" || x[i] == "Primary"){
y[i] <- "Until primary"
} else if (x[i] == "Secondary" || x[i] == "Tertiary"){
y[i] <- "Secondary or higher"
} else {
print(x[i]); y[i] <- x[i]
}
}
head(y, 10); table(y, useNA = "always")## [1] "Secondary or higher" "Secondary or higher" "Until primary" "Until primary" "Secondary or higher" "Secondary or higher" "Secondary or higher" "Until primary" "Secondary or higher" "Until primary"
## y
## Secondary or higher Unknown Until primary <NA>
## 1958 5 2790 0
z <- ifelse(x == "Secondary" | x == "Tertiary", "Secondary or higher", "Until primary")
head(z, 10); table(z, useNA = "always")## [1] "Secondary or higher" "Secondary or higher" "Until primary" "Until primary" "Secondary or higher" "Secondary or higher" "Secondary or higher" "Until primary" "Secondary or higher" "Until primary"
## z
## Secondary or higher Until primary <NA>
## 1958 2790 5
w <- case_when(is.na(x) ~ "Unknown",
x %in% c("Less primary", "Primary") ~ "Until primary",
x %in% c("Secondary", "Tertiary") ~ "Secondary or higher",
TRUE ~ x)
head(w, 10); table(w, useNA = "always")## [1] "Secondary or higher" "Secondary or higher" "Until primary" "Until primary" "Secondary or higher" "Secondary or higher" "Secondary or higher" "Until primary" "Secondary or higher" "Until primary"
## w
## Secondary or higher Unknown Until primary <NA>
## 1958 5 2790 0
## tibble [4,753 × 2] (S3: tbl_df/tbl/data.frame)
## $ score_lon1: num [1:4753] 3 3 3 3 3 3 6 3 3 6 ...
## $ score_sup1: num [1:4753] 12 12 12 11 12 13 10 12 11 13 ...
## $score_lon1
## [1] "numeric"
##
## $score_sup1
## [1] "numeric"
## score_lon1 score_sup1
## "numeric" "numeric"
## num [1:4753] 15 15 15 14 15 16 16 15 14 19 ...
## score_lon1 score_sup1
## 16712 50008
## [1] "1 4" "2 3" "3 2" "4 1"
* Del mismo modo que la función c, sapply intenta forzar a todos los elementos de la lista inicialmente resultante a tener el mismo tipo, en el orden NULL < lógico < entero < numérico < carácter (esencialmente)
## [1] "q0002_hhid" "sex" "age" "mar_1" "edu1" "phys_hea1" "hea1" "dep1" "score_lon1" "score_sup1" "income1" "income_inf1"
## [1] 4753 12
## [1] 12
## [1] 4753
## [1] "q0002_hhid" "dep2" "score_lon2" "score_sup2" "income2" "phys_hea2" "hea2"
## [1] 2400 7
## [1] "q0002_hhid" "dep3" "score_lon3" "score_sup3" "income3" "phys_hea3" "health3"
## [1] 1512 7
Cuántos missings tiene la variable q0002_hhid?
##
## FALSE
## 4753
##
## FALSE
## 2400
##
## FALSE
## 1512
Ya tenemos las bbdd abiertas y guardadas. Sólo tenemos que unir por la(s) variable(s) identificadora(s), en este caso q0002_hhid.
data <- dataw1 %>% # partimos de la base de datos dataw1 y entonces
full_join(dataw2, by = "q0002_hhid") %>% # unimos horizontalmente con todas las observaciones de dataw2 con q0002_hhid iguales a los de dataw1 y añadimos las nuevas, y entonces
full_join(dataw3, by = "q0002_hhid") # unimos horizontalmente con todas las observaciones de dataw3 con q0002_hhid iguales a los que ya había y añadimos las nuevas
names(data)## [1] "q0002_hhid" "sex" "age" "mar_1" "edu1" "phys_hea1" "hea1" "dep1" "score_lon1" "score_sup1" "income1" "income_inf1" "dep2" "score_lon2" "score_sup2" "income2"
## [17] "phys_hea2" "hea2" "dep3" "score_lon3" "score_sup3" "income3" "phys_hea3" "health3"
## [1] 4753 24
Abrimos las 3 olas de un ensayo clínico en fichero separados. Usamos la función read_dta en la librería haven.
encoding = "latin1" (información que proporciona con más detalle la ayuda de la función, ?read_dta)write_dta permite guardar data frames en formato de Stata para distintas versiones, especificando el argumento version.foreign y readstata13.library(haven)
ac1 <- read_dta("data_curs_stat/Assaig_clinic_w1.dta"); ac2 <- read_dta("data_curs_stat/Assaig_clinic_w2.dta"); ac3 <- read_dta("data_curs_stat/Assaig_clinic_w3.dta")Una base de datos en formato de Stata está guardada de manera distinta que una base de datos en formato de excel. Por ello, al abrirla como un data frame, éste tiene una estructura ligeramente distinta.
Recordemos que cuando abrimos la base de datos que guardamos en el objeto dataw1, la variable phys_hea1 era un vector carácter y sus elementos eran las distintas categorías, de manera que la transformamos de manera razonable en factor mediante la función as.factor. Ahora, sin embargo, al abrir bases de datos en formato Stata con la librería haven encontramos algunas diferencias
## [1] "haven_labelled" "vctrs_vctr" "double"
## dbl+lbl [1:400] 2, 1, 1, 3, 1, 2, 3, 2, 1, 1, 3, 2, 2, 1, 1, 3, 2, 1, 2, 2, 1, 1, 1, 1, 1, 1, 2, 3, 3, 3, 3, 1, 3, 1, 1, 2, 4, 3, 4, 3, 1, 2, 1, 2, 3, 1, 3, 1, 2, 2, 1, 1, 3, 1, 1, 3, 1, 1, 1, 2, 1, 2, 2, 1, 1, 2, 3, 2, 1, 1, 1, 1, 1,...
## @ label : chr "Level of education"
## @ format.stata: chr "%12.0g"
## @ labels : Named num [1:4] 1 2 3 4
## ..- attr(*, "names")= chr [1:4] "Less primary" "Primary" "Secondary" "Tertiary"
## $label
## [1] "Level of education"
##
## $format.stata
## [1] "%12.0g"
##
## $class
## [1] "haven_labelled" "vctrs_vctr" "double"
##
## $labels
## Less primary Primary Secondary Tertiary
## 1 2 3 4
##
## 1 2 3 4
## 179 112 89 19
En formatos de fichero distintos de excel o .csv, como pueden ser los de Stata o SPSS, las categorías de las variables categóricas pueden estar grabadas como etiquetas, mientras los valores son numéricos.
## [1] "factor"
## Factor w/ 4 levels "1","2","3","4": 2 1 1 3 1 2 3 2 1 1 ...
## $levels
## [1] "1" "2" "3" "4"
##
## $class
## [1] "factor"
##
## 1 2 3 4
## 179 112 89 19
En estos casos, la función as_factor de la librería haven es más adecuada que as.factor para recuperar las etiquetas, ya que, a diferencia de esta última, tiene en cuenta los atributos de la variable original.
## [1] "factor"
## Factor w/ 4 levels "Less primary",..: 2 1 1 3 1 2 3 2 1 1 ...
## - attr(*, "label")= chr "Level of education"
## $levels
## [1] "Less primary" "Primary" "Secondary" "Tertiary"
##
## $class
## [1] "factor"
##
## $label
## [1] "Level of education"
##
## Less primary Primary Secondary Tertiary
## 179 112 89 19
## [1] 400 9
## [1] "q0002_hhid" "number_id" "ID_ECS" "sex" "age" "mar1" "edu1" "hea1" "grups"
## [1] 400 6
## [1] "q0002_hhid" "number_id" "ID_ECS" "grups" "hea2" "dep2"
## [1] 400 6
## [1] "q0002_hhid" "number_id" "ID_ECS" "grups" "hea3" "dep3"
Para poder combinarlas verticalmente, los nombres de las columnas que queremos combinar tienen que ser iguales. Para ello usamos otra función de la librería dplyr: rename
acr1 <- ac1 %>% # partimos de ac1 y entonces
rename(hea = hea1) # renombramos hea1 como hea
acr2 <- ac2 %>% # partimos de ac2 y entonces
rename(hea = hea2, dep = dep2) # renombramos hea2 como hea y dep2 como dea
acr3 <- ac3 %>%
rename(hea = hea3, dep = dep3) # mismo proceso para ac3
#names(ac1); names(ac2); names(ac3)Finalmente combinamos las filas con otra función de la librería dplyr: bind_rows
Nota: La función rbind de R hace esencialmente lo mismo, pero necesita que las 3 bases de datos tengan exactamente las mismas columnas. En caso en que esto no ocurre, como el presente, da un error. Además, tanto en el caso de rbind como de bind_rows, conviene que las columnas por las que se combina (las de igual nombre) tengan la misma clase, pues en caso contrario pueden ocurrir errores o comportamiento extraños.
Nota: Las funciones cbind y bind_cols de R y dplyr respectivamente combinan por columnas. Se diferencian de un merge en que no hay una columna “común” por la que fusionar, sino que se añaden las columnas de los distintos objetos, tal como están ordenadas en cada uno de ellos. Además, todas las columnas tienen que tener el mismo número de elementos.
Fusionamos horizontalmente las 3 bbdd originales de los ensayos clínicos. En este caso no las podemos combinar porque tienen distinto número de filas:
ach <- ac1 %>%
full_join(ac2, by = c("q0002_hhid", "number_id", "ID_ECS", "grups")) %>%
full_join(ac3, by = c("q0002_hhid", "number_id", "ID_ECS", "grups"))
ach %>%
head()## # A tibble: 6 × 13
## q0002_hhid number_id ID_ECS sex age mar1 edu1 hea1 grups hea2 dep2 hea3 dep3
## <dbl> <chr> <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl> <dbl+lbl> <dbl> <dbl+lbl> <dbl> <dbl+lbl>
## 1 2306 724002306 2306 2 [fem] 3 [50-64] 0 [No] 2 [Primary] 45.9 1 [Intervention] NA NA NA NA
## 2 3363 724003363 3363 2 [fem] 4 [65-79] 0 [No] 1 [Less primary] 23.3 1 [Intervention] 24.6 0 [No] NA NA
## 3 3228 724003228 3228 2 [fem] 3 [50-64] 1 [Yes] 1 [Less primary] 49.1 1 [Intervention] 35.2 0 [No] NA NA
## 4 2371 724002371 2371 2 [fem] 1 [18-34] 0 [No] 3 [Secondary] 57.5 1 [Intervention] NA NA NA NA
## 5 2452 724002452 2452 2 [fem] 2 [35-49] 0 [No] 1 [Less primary] 41.3 1 [Intervention] 44.3 0 [No] 31.6 1 [Yes]
## 6 2750 724002750 2750 2 [fem] 5 [80+] 0 [No] 2 [Primary] 19.5 1 [Intervention] NA NA NA NA
## q0002_hhid number_id ID_ECS sex age mar1 edu1 hea1 grups q0002_hhid number_id ID_ECS grups hea2 dep2 q0002_hhid number_id ID_ECS grups hea3 dep3
## 1 2306 724002306 2306 2 3 0 2 45.88945 1 48 724000048 48 0 46.31060 0 48 724000048 48 0 47.35691 0
## 2 3363 724003363 3363 2 4 0 1 23.27672 1 61 724000061 61 0 53.24683 0 61 724000061 61 0 67.96230 0
## 3 3228 724003228 3228 2 3 1 1 49.08676 1 90 724000090 90 1 NA NA 90 724000090 90 1 NA NA
## 4 2371 724002371 2371 2 1 0 3 57.50416 1 128 724000128 128 0 60.88645 0 128 724000128 128 0 47.54218 0
## 5 2452 724002452 2452 2 2 0 1 41.27082 1 181 724000181 181 1 NA NA 181 724000181 181 1 NA NA
## 6 2750 724002750 2750 2 5 0 2 19.50435 1 191 724000191 191 0 NA NA 191 724000191 191 0 NA NA
Las funciones de las librerías de la familia tidyverse tratan los nombres de las variables de un data frame de una manera distinta a como lo hacen las funciones de R base. Por ejemplo, admiten determinados nombres especiales, pero sin embargo no admiten variables distintas con el mismo nombre, como sí ocurre en el último ejemplo.
Cuando se combinan las filas mediante bind_cols, las columnas que sólo están en una base de datos se llenan con missings. Por ejemplo, es el caso de las variables sociodemográficas, que en el dataframe acv sólo tienen datos en las filas correspondientes a la ola 1, que son las filas obtenidas de ac1.
## # A tibble: 6 × 11
## wave q0002_hhid number_id ID_ECS sex age mar1 edu1 hea grups dep
## <chr> <dbl> <chr> <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl> <dbl+lbl> <dbl+lbl>
## 1 1 48 724000048 48 2 [fem] 3 [50-64] 0 [No] 3 [Secondary] 47.9 0 [Control] NA
## 2 2 48 724000048 48 NA NA NA NA 46.3 0 [Control] 0 [No]
## 3 3 48 724000048 48 NA NA NA NA 47.4 0 [Control] 0 [No]
## 4 1 61 724000061 61 2 [fem] 2 [35-49] 0 [No] 4 [Tertiary] 63.2 0 [Control] NA
## 5 2 61 724000061 61 NA NA NA NA 53.2 0 [Control] 0 [No]
## 6 3 61 724000061 61 NA NA NA NA 68.0 0 [Control] 0 [No]
Acudimos ahora a la función fill de la librería tidyr:
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
##
## extract
acv %>%
group_by(q0002_hhid) %>% # agrupamos por id y entonces
arrange(wave) %>% # ordenamos por ola para tener la fila con datos antes que las otras, y entonces
fill(sex, age, mar1, edu1) %>% # para cada id rellenamos las filas vacías de las variables especificadas con los valores que no son missings, y entonces
ungroup() %>% # desagrupamos y entonces
arrange(q0002_hhid) %>% # ordenamos por id para ver las mismas filas que arriba, y entonces
head() # nos quedamos con las primeras filas## # A tibble: 6 × 11
## wave q0002_hhid number_id ID_ECS sex age mar1 edu1 hea grups dep
## <chr> <dbl> <chr> <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl> <dbl+lbl> <dbl+lbl>
## 1 1 48 724000048 48 2 [fem] 3 [50-64] 0 [No] 3 [Secondary] 47.9 0 [Control] NA
## 2 2 48 724000048 48 2 [fem] 3 [50-64] 0 [No] 3 [Secondary] 46.3 0 [Control] 0 [No]
## 3 3 48 724000048 48 2 [fem] 3 [50-64] 0 [No] 3 [Secondary] 47.4 0 [Control] 0 [No]
## 4 1 61 724000061 61 2 [fem] 2 [35-49] 0 [No] 4 [Tertiary] 63.2 0 [Control] NA
## 5 2 61 724000061 61 2 [fem] 2 [35-49] 0 [No] 4 [Tertiary] 53.2 0 [Control] 0 [No]
## 6 3 61 724000061 61 2 [fem] 2 [35-49] 0 [No] 4 [Tertiary] 68.0 0 [Control] 0 [No]
full_join creamos una base de datos resultado de fusionar las 3 iniciales e incluír todas las observaciones de cada una de ellas. Consultando en la ayuda, este ejercicio consiste en fusionar las 3 bases de datos, pero incluyendo sólo aquellas observaciones comunes a las 3. Cuántas observaciones tiene?acv; para qué se añadió .id = "wave"?; ejecutar rbind(ac1,ac2,ac3)cbind por bind_cols en el último ejemplo. Qué ocurre si se sustituye el valor por defecto del argumento .name_repair por los otros valores posibles?acv de tal manera que sus variables sociodemográficas estén completas tal como se indicó en la diapositiva anterior.achacv.Pivotar es el proceso de transformar una base de datos vertical/longitudinal a una horizontal o a la inversa.
Pivotaje usando la función pivot_wider de la librería tidyr:
#?pivot_wider
acv %>%
pivot_wider(names_from = c("wave"), values_from = c("hea", "dep")) %>%
# select(-where(~all(is.na(.)))) %>% # quitamos aquellas columnas que cumplen que todos sus valores son missings
head()## # A tibble: 6 × 14
## # Groups: q0002_hhid [6]
## q0002_hhid number_id ID_ECS sex age mar1 edu1 grups hea_1 hea_2 hea_3 dep_1 dep_2 dep_3
## <dbl> <chr> <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl> <dbl> <dbl> <dbl+lbl> <dbl+lbl> <dbl+lb>
## 1 2306 724002306 2306 2 [fem] 3 [50-64] 0 [No] 2 [Primary] 1 [Intervention] 45.9 NA NA NA NA NA
## 2 3363 724003363 3363 2 [fem] 4 [65-79] 0 [No] 1 [Less primary] 1 [Intervention] 23.3 24.6 NA NA 0 [No] NA
## 3 3228 724003228 3228 2 [fem] 3 [50-64] 1 [Yes] 1 [Less primary] 1 [Intervention] 49.1 35.2 NA NA 0 [No] NA
## 4 2371 724002371 2371 2 [fem] 1 [18-34] 0 [No] 3 [Secondary] 1 [Intervention] 57.5 NA NA NA NA NA
## 5 2452 724002452 2452 2 [fem] 2 [35-49] 0 [No] 1 [Less primary] 1 [Intervention] 41.3 44.3 31.6 NA 0 [No] 1 [Yes]
## 6 2750 724002750 2750 2 [fem] 5 [80+] 0 [No] 2 [Primary] 1 [Intervention] 19.5 NA NA NA NA NA
Si tenemos un sólo conjunto de columnas longitudinal, por ejemplo hea1, hea2 y hea3, podemos acudir a la función pivot_longer de la librería tidyr:
ach %>%
select(-starts_with("dep")) %>% # quitamos las columnas que empiezan con dep y entonces
pivot_longer(cols = starts_with("hea"), names_to = "wave", values_to = "hea") %>% # pivotamos, y entonces
arrange(q0002_hhid) %>% head()## # A tibble: 6 × 10
## q0002_hhid number_id ID_ECS sex age mar1 edu1 grups wave hea
## <dbl> <chr> <dbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <dbl+lbl> <chr> <dbl>
## 1 48 724000048 48 2 [fem] 3 [50-64] 0 [No] 3 [Secondary] 0 [Control] hea1 47.9
## 2 48 724000048 48 2 [fem] 3 [50-64] 0 [No] 3 [Secondary] 0 [Control] hea2 46.3
## 3 48 724000048 48 2 [fem] 3 [50-64] 0 [No] 3 [Secondary] 0 [Control] hea3 47.4
## 4 61 724000061 61 2 [fem] 2 [35-49] 0 [No] 4 [Tertiary] 0 [Control] hea1 63.2
## 5 61 724000061 61 2 [fem] 2 [35-49] 0 [No] 4 [Tertiary] 0 [Control] hea2 53.2
## 6 61 724000061 61 2 [fem] 2 [35-49] 0 [No] 4 [Tertiary] 0 [Control] hea3 68.0
Cuando hay más de un conjunto de columnas que queremos pivotar, como en este caso, en el que tenemos por un lado hea1, hea2 y hea3 y por otro lado, dep2 y dep3, usaremos la función reshape de R. Para ello, primero necesitamos el mismo número de columnas de hea y dep, por tanto creamos la columna dep1.
ach$dep1 <- NA # como en un data frame todas las variables tienen las mismas filas, R ya entiende que tiene que asignar NA para cada fila.También necesitamos que el objeto sea un data frame propiamente dicho, en lugar de un tibble (son así llamados los data frames creados por las funciones de la familia tidyverse, que se caracterizan porque su clase es c("tbl_df", "tbl", "data.frame"))
## [1] "tbl_df" "tbl" "data.frame"
## [1] "data.frame"
reshape(as.data.frame(ach), varying = list(paste0("hea", 1:3), paste0("dep", 1:3)), v.names = c("hea", "dep"), timevar = "wave", idvar = "q0002_hhid", direction = "long") %>%
arrange(q0002_hhid) %>%
head()## q0002_hhid number_id ID_ECS sex age mar1 edu1 grups wave hea dep
## 48.1 48 724000048 48 2 3 0 3 0 1 47.94347 NA
## 48.2 48 724000048 48 2 3 0 3 0 2 46.31060 0
## 48.3 48 724000048 48 2 3 0 3 0 3 47.35691 0
## 61.1 61 724000061 61 2 2 0 4 0 1 63.24433 NA
## 61.2 61 724000061 61 2 2 0 4 0 2 53.24683 0
## 61.3 61 724000061 61 2 2 0 4 0 3 67.96230 0
##
## Pearson's product-moment correlation
##
## data: dataw1$score_sup1 and dataw1$income1
## t = 3.2915, df = 3518, p-value = 0.001007
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.02241180 0.08828366
## sample estimates:
## cor
## 0.05540802
##
## Shapiro-Wilk normality test
##
## data: dataw1$hea1
## W = 0.98634, p-value < 2.2e-16
## Warning in ks.test(dataw1$hea1, "pnorm"): ties should not be present for the Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: dataw1$hea1
## D = 0.99954, p-value < 2.2e-16
## alternative hypothesis: two-sided
Ambos test son significativos, así pues descartamos la hipótesis nula consistente en la normalidad de la distribución de hea1 para la muestra en dataw1.
Para otros tests de normalidad como los de Lilliefors, Anderson-Darling o Jarque-Bera, disponemos de funciones en librerías como DescTools, nortest o tseries.
##
## fem masc
## No 1282 626
## Yes 1320 1525
test <- chisq.test(dataw1$mar_1, dataw1$sex)
test # por defecto, test con correccion de continuidad si es 2x2##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: dataw1$mar_1 and dataw1$sex
## X-squared = 198.48, df = 1, p-value < 2.2e-16
## dataw1$sex
## dataw1$mar_1 fem masc
## No 1044.523 863.4774
## Yes 1557.477 1287.5226
##
## Pearson's Chi-squared test
##
## data: dataw1$mar_1 and dataw1$sex
## X-squared = 199.31, df = 1, p-value < 2.2e-16
##
## Fisher's Exact Test for Count Data
##
## data: dataw1$mar_1 and dataw1$sex
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 2.093293 2.674704
## sample estimates:
## odds ratio
## 2.365541
## [1] 0.00137449
##
## One Sample t-test
##
## data: dataw1$new_var
## t = 0.092092, df = 4752, p-value = 0.9266
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.02788569 0.03063467
## sample estimates:
## mean of x
## 0.00137449
t.test(dataw1$new_var)$conf.int # intervalo de confianza para la media muestral de new_var en dataw1## [1] -0.02788569 0.03063467
## attr(,"conf.level")
## [1] 0.95
## [1] 53.96553
##
## Wilcoxon signed rank test with continuity correction
##
## data: dataw1$hea1
## V = 5409216, p-value = 0.0278
## alternative hypothesis: true location is not equal to 54
## 95 percent confidence interval:
## 54.06244 55.07854
## sample estimates:
## (pseudo)median
## 54.57055
Como conclusiones, no podemos descartar que la media de new_var sea 0, pero sí podemos descartar que los valores de hea1 se localizen en torno a 54.
Comparamos entre grupos mediante un t-test:
##
## Welch Two Sample t-test
##
## data: hea1 by sex
## t = -14.557, df = 4523.6, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group fem and group masc is not equal to 0
## 95 percent confidence interval:
## -7.861186 -5.995119
## sample estimates:
## mean in group fem mean in group masc
## 50.82465 57.75280
t.test(data$hea1, data$hea2, paired = TRUE) # comparamos el estado de salud en la ola 1 con la ola 2##
## Paired t-test
##
## data: data$hea1 and data$hea2
## t = 0.81332, df = 2384, p-value = 0.4161
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2407952 0.5820965
## sample estimates:
## mean of the differences
## 0.1706506
t.test(hea1 ~ sex, data = dataw1, mu = -7) # comparamos la diferencia de medias con el valor mu (hipótesis nula)##
## Welch Two Sample t-test
##
## data: hea1 by sex
## t = 0.15097, df = 4523.6, p-value = 0.88
## alternative hypothesis: true difference in means between group fem and group masc is not equal to -7
## 95 percent confidence interval:
## -7.861186 -5.995119
## sample estimates:
## mean in group fem mean in group masc
## 50.82465 57.75280
Nota: mediante el argumento alternative de la función t.test, podemos especificar si la hipótesis nula es la igualdad (test bilateral) o una desigualdad.
wilcox.test(hea1 ~ sex, data = dataw1,
alternative = "two.sided", mu = 0,
paired = FALSE, exact = FALSE,
conf.int = TRUE, conf.level = 0.95)##
## Wilcoxon rank sum test with continuity correction
##
## data: hea1 by sex
## W = 1964188, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
## -7.958215 -6.019921
## sample estimates:
## difference in location
## -6.985333
Descartamos la hipótesis nula, por tanto podemos concluír que la distribución de hea1 en los hombres es significativamente diferente de la distribución en las mujeres.
##
## One-sample Kolmogorov-Smirnov test
##
## data: dataw1$new_var
## D = 0.013135, p-value = 0.3851
## alternative hypothesis: two-sided
## Df Sum Sq Mean Sq F value Pr(>F)
## edu1 3 2 0.631 0.596 0.618
## Residuals 4744 5026 1.059
## 5 observations deleted due to missingness
##
## Call:
## lm(formula = hea1 ~ age + sex + income1 + dep1, data = dataw1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -53.560 -8.844 0.710 9.444 40.641
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 64.99142 0.78354 82.946 < 2e-16 ***
## age35-49 -5.41946 0.95045 -5.702 1.28e-08 ***
## age50-64 -11.26679 0.79368 -14.196 < 2e-16 ***
## age65-79 -19.19973 0.80852 -23.747 < 2e-16 ***
## age80+ -30.15751 1.02706 -29.363 < 2e-16 ***
## sexmasc 5.07076 0.44117 11.494 < 2e-16 ***
## income1 0.16659 0.02742 6.075 1.37e-09 ***
## dep1Yes -14.52514 0.67309 -21.580 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.09 on 3667 degrees of freedom
## (1078 observations deleted due to missingness)
## Multiple R-squared: 0.3837, Adjusted R-squared: 0.3825
## F-statistic: 326.1 on 7 and 3667 DF, p-value: < 2.2e-16
## (Intercept) age35-49 age50-64 age65-79 age80+ sexmasc income1 dep1Yes
## 64.991424 -5.419455 -11.266791 -19.199732 -30.157514 5.070758 0.166592 -14.525141
## 2.5 % 97.5 %
## (Intercept) 63.4552039 66.5276443
## age35-49 -7.2829188 -3.5559915
## age50-64 -12.8228922 -9.7106894
## age65-79 -20.7849349 -17.6145292
## age80+ -32.1711810 -28.1438463
## sexmasc 4.2058013 5.9357156
## income1 0.1128245 0.2203594
## dep1Yes -15.8448166 -13.2054651
#?glm
fit <- glm(dep1 ~ age + sex + income1*hea1 + mar_1*score_lon1, data = dataw1, family = binomial)
summary(fit)##
## Call:
## glm(formula = dep1 ~ age + sex + income1 * hea1 + mar_1 * score_lon1,
## family = binomial, data = dataw1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0231 -0.4485 -0.2799 -0.1678 3.2085
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9027594 0.4760209 1.896 0.057898 .
## age35-49 0.1435643 0.3101981 0.463 0.643497
## age50-64 -0.2696730 0.2772296 -0.973 0.330681
## age65-79 -1.0277636 0.2922478 -3.517 0.000437 ***
## age80+ -2.0710935 0.3598878 -5.755 8.67e-09 ***
## sexmasc -0.2519979 0.1314381 -1.917 0.055208 .
## income1 -0.0154454 0.0276949 -0.558 0.577051
## hea1 -0.0720689 0.0064107 -11.242 < 2e-16 ***
## mar_1Yes -0.5831916 0.3266439 -1.785 0.074196 .
## score_lon1 0.3268351 0.0406748 8.035 9.33e-16 ***
## income1:hea1 -0.0001597 0.0005689 -0.281 0.778968
## mar_1Yes:score_lon1 0.1290971 0.0668269 1.932 0.053383 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2585.4 on 3566 degrees of freedom
## Residual deviance: 1937.0 on 3555 degrees of freedom
## (1186 observations deleted due to missingness)
## AIC: 1961
##
## Number of Fisher Scoring iterations: 6
## (Intercept) age35-49 age50-64 age65-79 age80+ sexmasc income1 hea1 mar_1Yes score_lon1 income1:hea1 mar_1Yes:score_lon1
## 2.4663996 1.1543810 0.7636291 0.3578063 0.1260479 0.7772463 0.9846733 0.9304668 0.5581142 1.3865729 0.9998403 1.1378006
ggplot: una única variable predictorafit <- lm(hea1 ~ income1, data = dataw1)
dataw1 %>%
filter(!is.na(hea1), !is.na(income1)) %>% # nos quedamos con las observaciones que no tienen missings en el modelo
# filter(!if_any(c(hea1, income1), is.na)) %>% # lo mismo de forma resumida
cbind(predict(fit, interval = "confidence")) %>% # añadimos la predicción del modelo y los intervalos de confianza como nuevas columnas
ggplot(aes(x = income1)) + # creamos el cuadro y añadimos el eje x
geom_point(aes(y = hea1)) + # añadimos los punos d la base de datos
geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = .15) + # añadimos los CI; también es posible como sigue
# geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
# geom_line(aes(y = upr), color = "red", linetype = "dashed") +
annotate("text",
label = paste("La correlación es", round(cor(dataw1$hea1, dataw1$income1, use = "complete.obs"), 2)),
x = 30, y = 15) +
theme_bw() # cambiamos el tema por defectoggplot: una variable predictora continua y otra categóricalibrary(ggplot2)
fit <- lm(hea1 ~ income1 + edu1, data = dataw1)
dataw1 %>%
# filter(!is.na(hea1), !is.na(income1), !is.na(edu1)) %>% # nos quedamos con las observaciones que no tienen missings en el modelo
filter(!if_any(c(hea1, income1, edu1), is.na)) %>% # lo mismo de forma resumida
cbind(predict(fit, interval = "confidence")) %>% # añadimos la predicción del modelo y los intervalos de confianza como nuevas columnas
ggplot(aes(x = income1, color = edu1)) + # creamos el cuadro y añadimos el eje x
geom_point(aes(y = hea1)) + # añadimos los punos d la base de datos
geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
geom_ribbon(aes(ymin = lwr, ymax = upr, color = NULL, fill = edu1), alpha = .15) # añadimos los CIggplot: una única variable predictoralibrary(ggplot2)
fit <- glm(dep1 ~ income1, data = dataw1, family = binomial)
dataw1 %>%
filter(!if_any(c(dep1, income1), is.na)) %>% # nos quedamos con las observaciones que no tienen missings en el modelo
mutate(fit = predict(fit, type = "response")) %>% # añadimos la predicción del modelo como nuevas columnas
ggplot(aes(x = income1)) + # creamos el cuadro y añadimos el eje x
geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
labs(y = "Predicted probability of dep1 == Yes",
title = "Predicting depression from income through a logistic model",
caption = "Data source: A cohort study") + # especificamos la etiqueta del eje y, del título y de pie de figura
theme_classic() #cambiamos el tema por defectoggplot: una variable predictora continua y otra categóricalibrary(ggplot2)
fit <- glm(dep1 ~ income1 + edu1, data = dataw1, family = binomial)
dataw1 %>%
filter(!if_any(c(dep1, income1, edu1), is.na)) %>% # nos quedamos con las observaciones que no tienen missings en el modelo
mutate(fit = predict(fit, type = "response")) %>% # añadimos la predicción del modelo como nuevas columnas
ggplot(aes(x = income1, color = edu1)) + # creamos el cuadro y añadimos el eje x
geom_line(aes(y = fit)) + # añadimos la línea que el modelo induce
labs(y = "Predicted probability of dep1 == Yes") + # especificamos la etiqueta del eje y
theme_minimal() # cambiamos el tema por defectoDe la misma manera que tidyverse hace referencia a una colección de librerías que simplifican la manipulación de datframes, existe otra colección de librerías, easystats (https://easystats.github.io/easystats/), cuyo objetivo es mostrar los resultados de tests y modelos estadísticos con mejor formato.
La librería modelsummary dispone de 2 funciones que permiten visualizar multiples resultados en tablas bien organizadas:
modelsummary (https://vincentarelbundock.github.io/modelsummary/articles/modelsummary.html) ydatasummary (https://vincentarelbundock.github.io/modelsummary/articles/datasummary.html).Otras librerías que permiten visualizar descriptivos y resultados resumidos en tablas son:
summarytools (https://cran.r-project.org/web/packages/summarytools/vignettes/introduction.html) ycompareGroups (https://cran.r-project.org/web/packages/compareGroups/vignettes/compareGroups_vignette.html).Existen distintas librerías que facilitan la representación gráfica de modelos, predicciones, …:
jtools (https://cran.r-project.org/web/packages/jtools/vignettes/effect_plot.html)effects (https://cran.r-project.org/web/packages/effects/vignettes/predictor-effects-gallery.pdf)ggeffects(https://strengejacke.github.io/ggeffects/)Ejercicios:
cor. Calcular la correlación de Spearman. Cómo tendríamos que hacer si hay valores missing?performance y parameters y probar las funciones model_performance y model_parameters aplicadas a fit.## [1] 4
## [1] "x"
## [1] 4
## [1] 5
.r), se puede ejecutar a través de la función source. Por ejemplo source("fichero.r").nchar.## [1] TRUE
haven guarda las etiquetas de los objetos como atributos de los mismos, a los que se puede acceder por medio de la función attributes otras librerías como labelled, sjmisc o sjlabelled también tienen múltiples funciones creadas para trabajar con etiquetas.sample permite seleccionar elementos de un vector aleatoriamente. La función slice_sample de la librería dplyr permite seleccionar filas (observaciones) al azar de un data frame.ggsave permite guardar gráficas generadas a partir de ggplot, especificando el formato, medida y calidad entre otros.plyr y dplyr comparten varias funciones, entre ellas rename. Si se carga primero dplyr y después plyr R asocia rename a la función de plyr, pero al usarla como la función de dplyr, es probable que se produzca un error.::##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise, summarize
# ?rename # con la ayuda podemos ver que esta función está en varias librerías actualmente cargadas
names(ac1)## [1] "q0002_hhid" "number_id" "ID_ECS" "sex" "age" "mar1" "edu1" "hea1" "grups"
## Error in rename(ac1, hea = hea1): unused argument (hea = hea1)
## [1] "q0002_hhid" "number_id" "ID_ECS" "sex" "age" "mar1" "edu1" "hea" "grups"
names(plyr::rename(ac1, hea = hea1)) # comprobamos que la función que estaba usando antes era la de plyr, ya que genera el mismo error## Error in plyr::rename(ac1, hea = hea1): unused argument (hea = hea1)
rename <- dplyr::rename # otra solución es asignar al nombre la función de dplyr
names(rename(ac1, hea = hea1)) # ya no hay error## [1] "q0002_hhid" "number_id" "ID_ECS" "sex" "age" "mar1" "edu1" "hea" "grups"
Tutoriales introductorios a R:
Blogs:
psych):
Multitud de libros:
Cursos online (MOOC’s) en plataformas como Coursera, edX, etc.
Acerca de librerías o materias específicas:
ggplot2ggplot2Recursos de RStudio
tidyverse o de actualizaciones relevantesPodéis consultarme vustras dudas a través de mi dirección de email: iago.gine@sjd.es